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ABSTRACT 


EP] 10 Jan 2022 


. Context. Ultra-short period planets undergo strong tidal interactions with their host star which lead to planet deformation and orbital tidal decay. 
-Cd Aims. WASP-103b is the exoplanet with the highest expected deformation signature in its transit light curve and one of the shortest expected 
Cospiral-in times. Measuring the tidal deformation of the planet would allow us to estimate the second degree fluid Love number and gain insight 
! into the planet's internal structure. Moreover, measuring the tidal decay timescale would allow us to estimate the stellar tidal quality factor, which 
is key to constraining stellar physics. 
Methods. We obtained 12 transit light curves of WASP-103b with the CHaracterising ExOplanet Satellite (CHEOPS) to estimate the tidal defor- 
mation and tidal decay of this extreme system. We modelled the high-precision CHEOPS transit light curves together with systematic instrumental 
„` 7 noise using multi-dimensional Gaussian process regression informed by a set of instrumental parameters. To model the tidal deformation, we used 
a parametrisation model which allowed us to determine the second degree fluid Love number of the planet. We combined our light curves with 
— previously observed transits of WASP-103b with the Hubble Space Telescope (HST) and Spitzer to increase the signal-to-noise of the light curve 
and better distinguish the minute signal expected from the planetary deformation. 
Results. We estimate the radial Love number of WASP-103b to be hy = 1.597055. This is the first time that the tidal deformation is directly 
(N detected (at 3 0) from the transit light curve of an exoplanet. Combining the transit times derived from CHEOPS, HST, and Spitzer light curves 
with the other transit times available in the literature, we find no significant orbital period variation for WASP-103b. However, the data show a hint 
of an orbital period increase instead of a decrease, as is expected for tidal decay. This could be either due to a visual companion star if this star is 
bound, the Applegate effect, or a statistical artefact. 

* Conclusions. The estimated Love number of WASP-103b is similar to Jupiter's. This will allow us to constrain the internal structure and compo- 
™= sition of WASP-103b, which could provide clues on the inflation of hot Jupiters. Future observations with James Webb Space Telescope (JWST) 
can better constrain the radial Love number of WASP-103b due to their high signal-to-noise and the smaller signature of limb darkening in the 
infrared. A longer time baseline is needed to constrain the tidal decay in this system. 
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- 
© 1. Introduction et al.[2021) Guaranteed Time Observing (GTO) programme, we 


are investigating the tidal interaction between ultra-hot Jupiters 


The extreme environment that ultra-short orbital period planets and their parent stars by attempting to measure their tidal decay 
are subjected to makes them ideal laboratories to study plan- — and deformation. 


etary physics. In addition to the very high temperatures, they 

also suffer from intense tidal forces which lead to a deformation Tidal forces tend to circularise planetary orbits and to syn- 
of the planet's shape and shrink-  chronise the planetary and stellar rotation with the orbital period. 
age of the planet's orbit. Hence, their study allows us to gain In hot Jupiter systems, the orbits are usually circularised and the 
a wealth of information on planet-to-star tidal interactions. As planet rotation is synchronised (Ogilvie & Lin|2004). However, 
part of the CHaracterising ExOplanet Satellite (CHEOPS) the synchronisation of the stellar rotation 1s still incomplete due 
to the longer and still poorly unconstrained timescale of this pro- 


* E-mail: susana.barros @astro.up.pt cess. (1980) showed that for planets with an orbital period 
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shorter than a third of the rotation period of the star, as it turns 
out to be the case for hot Jupiters, the tidal interaction leads to the 
unstable transfer of angular momentum from the planetary orbit 
to the stellar angular momentum. This results in the planet spi- 
ralling inwards and eventually being engulfed by the star. There- 
fore, tidal interactions between a star and a close-in exoplanet 
lead to shrinkage of the orbit and eventual tidal disruption of the 
planet. The synchronisation timescale of the stellar rotation de- 
pends on the tidal quality factor Q’, which is poorly constrained. 

The parameter Q' allows us to constrain stellar physics (e.g. 
Ogilvie & Lin/2007) and hence many attempts have been made 
to measure it. Studies of binary stars have estimated the tidal 
quality factor to be between 10° — 107 
(2005). However, hot Jupiter systems could be in a different tidal 
regime with a higher tidal factor (Q^ = 10°) and a weaker tidal 
decay (Ogilvie & Linp007| Penev & SasselovD01T). Estimates 
of Q’, through the measurement of the orbital period decrease 
were successful for the following hot Jupiters: WASP-4 (Q' = 


10* 2019) and WASP-12 (Q; = 10? 
2016$ |Yee et al.|2020). However, the measured values of 


Q; are lower than expected by theory (implying a stronger tidal 
dissipation) and it has not been possible to completely rule out 
other causes for the period decrease in these systems, such as ap- 
sidal precession. Statistical studies of the ensemble of known hot 
Jupiters show two regimes of tidal dissipation strength. The ma- 
jority of the studied systems had logio Q, = 8.26 + 0.14, while 


a smaller group had logio Q; = 7.3 + 0.4 (Collier Cameron &| 
Jardine018). 

The tidal deformation of a planet mostly depends on the 
planet-to-star distance and it is most significant for large planets 
that are almost filling their Roche lobe (e.g. [Ferraz-Mello et al.| 
[2008). Hence, it is larger for ultra-hot Jupiters. The radial defor- 
mation of a planet due to a perturbing potential can be quanti- 
fied using the second degree fluid Love number hy Es EAN, 
The Love number measures the distribution of mass within the 
planet depending on the concentration of heavy elements in the 
core of the planet relative to the envelope of the planet. There- 
fore, it provides insight into the internal structure differentiation 
of planets (Kramm et al.|2011). (2014) shows that hy 
is proportional to an asymmetry parameter q which relates the 
three axes (ri, r2, and r3) of the planetary ellipsoidal shape — if 
ry = ro(1 + 3q) and r5 = r;(1 — q), then 


a) 


with the volumetric radius Ry = ¥rir2r3 M, and M, being the 
planetary and the stellar mass, respectively, R, being the stel- 
lar radius, and a being the semi-major axis of the planet's or- 
bit. [Correia] (2014) also shows that the non-spherical shape of a 
deformed planet along with its varying projected area during a 
transit modifies the transit light curve and causes anomalies in 
the ingress, egress, and mid-transit phases compared to a spheri- 


cal case (Figure A1 of|Correia|2014). Detecting the deformation- 


induced signature in the light curve can therefore allow for the 
measurement of the Love number 
lard et al.[2019). 

Of the several attempts made to measure the deformation sig- 
nature, the most constraining is for WASP-121 using two Hubble 
Space Telescope (HST) /Space Telescope Imaging Spectrograph 
(STIS) transits (hy = 1.39 + 0.8 — < 2c significance d 
[lard et al.[2019). A measurement of an exoplanet's Love number 
was made for HAT-P-13b (Buhler et al.[2016). HAT-P-13b has a 
unique orbital configuration that allows for the measurement of 
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the Love number using apsidal precession. However, in this case 
some assumptions are required to estimate the Love number (see 
Section [5.4}. constrained the Love number 
of HAT-P-13b to be 1.116 < hy < 1.425. This was later up- 
dated to hy = 1.3 aie and allowed constraints on the maximum 
core size and the metallicity of the planet’s envelope, showing 
the power of the Love number in providing insights into the in- 
ternal structure of the planet 
(2016). Furthermore, unveiling the internal structure is in turn im- 
portant for understanding the formation of the planet itself since 
the distribution of the heavy elements and the core mass directly 
depend on formation mechanisms (Mordasini et al.[2012). 

Taking advantage of the high-precision and high pointing 
flexibility of the CHEOPS satellite, we designed a programme 
to measure the tidal decay and deformation of ultra-hot Jupiters. 
The expected amplitude of the deformation signature is largest 
for WASP-103b (~ 60 ppm) due to its larger radius among the 
ultra-hot Jupiters. Hence, this target was a priority for our pro- 
gramme. WASP-103b is a 1.5 Myyp and 1.5 Ryyp planet in a 22 
hour orbit around a late F-type star with a G magnitude of 12.2 
(Gillon et al./2014). The small amplitude of the tidal deforma- 
tion signal has prevented its detection until now and requires that 
CHEOPS transits are combined with other high signal-to-noise 
transits in order to allow us to estimate the planet's Love num- 
ber. The required long baseline of observations to measure the 
tidal decay of exoplanets also requires that the derived transit 
times from CHEOPS are combined with previously derived tran- 
sit times. 

In this paper, we present the first results of our tidal decay 
and deformation programme targeting WASP- 103b. In Section 2 
we describe the CHEOPS observations and in Section 3 we de- 
scribe complementary observations necessary to better constrain 
the system. In Section 4, we present our results for the variation 
of the planetary orbital period and discuss possible scenarios to 
explain it. In Section 5, we present our modelling of the tidal 
deformation combining CHEOPS results with HST and Spitzer 
observations. Finally, we draw our conclusions in Section 6. 


2. Observations, data reduction, and analysis 
2.1. CHEOPS observations of WASP-103b 


The objective of CHEOPS is to achieve a detailed characteri- 
sation of known exoplanets through high-precision photomet- 
ric observations. It is the first S-class ESA mission and it was 
launched on 18 December 2019 with science 
observations starting in April 2020. We obtained data as part of 
the CHEOPS Guaranteed Time Observing (GTO) programme: 
"Tidal decay and deformation (ID 0013)’. This programme aims 
to measure the tidal deformation and decay of short period exo- 
planets in order to constrain the planetary Love number and the 
stellar tidal dissipation parameter. This programme is included in 
one of the six GTO themes called feature characterisation which 
also includes one programme to search for moons and rings and 
one programme to measure the angle between the planetary orbit 
and the stellar spin through the gravity darkening effect. 

Currently the tidal deformation programme includes the tar- 
gets WASP-12b and WASP-103b. These, together with WASP- 
121b, are the best known targets to measure the tidal deforma- 
tion directly from the light curve (Akinsanmi et al.[2019). Un- 
fortunately, WASP-121b is not observable by CHEOPS due to 
pointing restrictions. 

Due to the extremely high photometric precision necessary 
to measure the tidal deformation, the original plan was to obtain 
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20 transits per year over the 3 years of the GTO. Due to the best 
visibility and observational efficiency of WASP-103 compared 
with WASP-12, this target was given priority and 20 transits 
were requested in the first year of the CHEOPS nominal mission. 
CHEOPS data suffer from interruptions due to Earth occultations 
or passages through the South Atlantic Anomaly (SAA), which 
can affect our observations. Therefore, extensive tests were per- 
formed during the preparation of the GTO of CHEOPS that show 
that a good coverage of ingress and egress is crucial in order to 
obtain accurate and precise times and also to better sample the 
shape deformation. 

We requested 90% efficiency in ingress and egress and 60% 
overall efficiency of transit observations. Since this would de- 
crease the number of possible observable transits, we started by 
requesting 90% efficiency in either ingress or egress. However, 
the first three transits showed poorer precision of the derived 
transit times, and hence, we included a stronger constraint of 
having 90% efficiency in both ingress and egress. This allowed 
us to observe only 12 of the 20 requested transits, but with in- 
creased accuracy for the derived transit times. After the first three 
test observations, we also increased the total requested time per 
observation. Originally, we requested the observations to cover 
three transit durations. For WASP-103, this corresponds to ~ 3.4 
CHEOPS orbits (~ 7.8 hours). However, for observations with 
an efficiency of less than 88%, we do not have the recommended 
three CHEOPS orbits of data to be able to detrend the system- 
atic noise (Maxted et al.|2021). Hence, we increased the duration 
of the observations which resulted in a much better detrending 
of the systematic noise (see Sectionp.2). The observation log is 
presented in Table[T] 


2.2. Photometric extraction 


The CHEOPS observations were reduced with the CHEOPS data 
reduction pipeline (DRP) (Hoyer et al.|2020). The DRP automat- 
ically processes all the CHEOPS data. It makes bias, dark, and 
flat corrections, and it applies gain, scattered light, and a cor- 
rection for the non-linearity of the detector response. The DRP 
simulates the field of view using the magnitudes and positions 
of stars in the Gaia DR2 catalogue 
[2018). These simulations are used to calculate the contamination 
of the target aperture by nearby stars. Due to the irregular PSF 
shape coupled with the rotation of the field of view of CHEOPS, 
the target star suffers from variable contamination from nearby 
stars. This contamination is a function of the angle of rotation of 
the satellite (roll angle) and the pointing jitter. The DRP calcu- 
lates and provides the contamination of the target aperture as a 
function of time so it can be corrected later. In the case of WASP- 
103, the simulation of the field of view shows a contaminating 
star inside the aperture ~ 16 arcsec from the target, as is shown 
in Figure [I] This contaminant adds ~ 0.9 % to the total flux in 
the aperture. 

The DRP also corrects the smearing trails of bright stars in 
the field of view. Due to the rotation of the field of view, this 
leads to a variable contamination of the target aperture. The DRP 
extracts the photometry for four apertures, three with a fixed ra- 
dius (22.5, 25, and 30 pixels) and one with an optimal radius, 
labelled RINF, DEFAULT, RSUP, and OPTIMAL. The radius of 
the optimal aperture is calculated for each data set to maximise 
the signal-to-noise ratio of the light curve. The DRP also cor- 
rects the background light which is estimated from an annulus 
around the target. The DRP produces a fits file with four ex- 
tracted light curves together with auxiliary information includ- 
ing, for example, the time series of the roll angle, the estimated 
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Fig. 1. Top panel: Field of view of WASP-103 as observed by CHEOPS 
transit #10. Bottom panel: Field of view of WASP-103 simulated by 
the DRP based on the Gaia catalogue and excluding WASP-103. The 
red cross marks the target position, while the red circle shows the DE- 
FAULT aperture. The image scale is 1 arcsec per pixel. 


contamination, the subtracted background, a quality flag, and the 
centroid position of the target star. These can be used to correct 
any systematic effects in the light curves. Furthermore, the DRP 
produces a report that states the performance of each step of the 
pipeline. More information about the data reduction pipeline can 
be found in [Hoyer et al.| (2020). For WASP-103, we considered 
both the OPTIMAL and DEFAULT apertures and chose the one 
with lower residuals in the final analysis as explained in the next 
sub-section. 


2.3. CHEOPS data analysis 


The light curves obtained with the DRP were corrected for the 
estimated contamination of the aperture in each light curve. 
Some of our observations were also affected by the atmospheric 
air glow at the beginning and end of an Earth occultation. In 
these cases, the air glow contaminates the observation, increas- 
ing the background to values higher than the target and ulti- 
mately leading to saturation. These points cannot be corrected 
for and we removed all of those with a background noise higher 
than the median noise of the target star. We also removed out- 
liers using 5 c clipping. In our case, all the light curves show a 
strong correlation with the roll angle (Figure [2). Moreover, the 
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Table 1. Log of CHEOPS observations of WASP-103b. 


# Start date Duration Noy;  Effic. APER Rap Decorrelation 
(UTC) (hours) (%) (pixels) 

1 | 2020-04-18T22:55:40.965587 6.02 269 74 OPTIMAL 17.0 roll , bg 

2  2020-05-02T19:53:40.996584 6.14 296 80 DEFAULT 25 roll 

3  2020-05-05T14:10:00.434364 6.55 308 78 OPTIMAL 17.5 bg 

4  2020-05-16T14:57:01.026698 9.64 523 90 DEFAULT 25 roll, x , y 

5  2020-05-19T10:08:00.408947 9.67 544 94 OPTIMAL 19.0 roll , bg 

6 | 2020-05-25T21:19:00.404138 9.97 560 94 OPTIMAL 19.0 roll, x 

7  2020-06-06T22:07:00.904021 9.37 540 96 DEFAULT 25.0 roll 

8  2020-06-07T20:04:00.500261 9.64 546 94 DEFAULT 25.0 roll 

9  2020-06-14T08:08:00.511905 9.64 533 92 OPTIMAL 19.0 roll , bg 

10 2020-06-18T23:06:00.996556 9.64 555 96 DEFAULT 25.0 roll 

11  2020-06-19T21:02:00.713493 9.64 538 93 OPTIMAL 19.0 roll , bg 

12  2020-06-20T19:07:39.419756 9.55 537 93 OPTIMAL 18.5 roll , cont, x, bg 


Notes. 


The transits are labelled by their sequence number throughout the paper. Effic. is the proportion of the time in which unobstructed observations of 
the target occurred. R,, is the aperture radius used for the photometric extraction. We also give the decorrelation parameters used for each light 
curve roll angle (roll), background (bg), x centroid (x), y centroid (y), and contamination (cont). 


light curves show extra correlations with a mix of instrumen- 
tal parameters, for example, the background flux, contamination 
rate, and x position. These light curves are presented in Figure[3] 
They clearly show systematic noise which is the residual of the 
variable contamination of the aperture, mentioned above, and it 
is highly correlated with instrumental parameters such as the roll 
angle of the satellite, the position of the star, and the background 
flux. 

During the preparation of the CHEOPS mission, several 
methods were tested to correct the systematic noise due to the 
rotating field of view and it was concluded that a better accuracy 
is achieved if the systematics are corrected simultaneously with 
the transit modelling. In order to derive transit parameters and si- 
multaneously decorrelate the CHEOPS light curves, we used two 
methods. The first one is based on multi-dimensional Gaussian 
processes (GPs) coupled with the 
batman transit model (Kreidberg|2015a). The second method is 
based on linear decorrelation using a combination of sinusoidal 


functions implemented in pycheops (Maxted et al.|2021). 


2.3.1. Multi-dimentional GP 


We performed GP regression with a Matern-3/2 kernel to model 
the flux dependence on the instrumental parameters using the 


George package (Ambikasaran et al.|2015). This is coupled 


with a transit model using a parametrisation and method simi- 
lar to the one used in [Barros et al.) (2020). The transit model is 
parameterised by the orbital period (P), mid-transit time (To), 
normalised separation of the planet (a/R,), the planet-to-star ra- 
dius ratio (r,/R,.), inclination (inc), and quadratic limb darken- 
ing law. We assumed a circular orbit 
[et al.[2018] and Section |3.1.2). The hyper-parameters of the GP 


are an amplitude (amp) and a length scale (s). 


For the shape parameters of the transit (a/R,, rp/Rx, and 
inc), we used Gaussian priors based on the results of 


(2018). For the limb darkening parameters, we also used 


Gaussian priors whose values and the uncertainties were de- 


rived with the LDTK code (Parviainen & Aigrain|2015 
2013) in the CHEOPS bandpass. For the mid-transit time, 


we used a uniform prior and we assumed the ephemerides de- 


rived by (2020) (To = 2456836.29630(07) and 
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Fig. 2. Residuals of the transit fit as a function of the roll angle of the 
satellite for transit number 8. A clear dependence is seen which is well 
fitted with a non-parametric GP model overplotted in orange (see Sec- 
tion [2.3.1]. The best fit length scale hyper-parameter in this case was 
TER degrees. The flux dependence with a roll angle is clearly seen in 
all our light curves, but the shape of the dependence varies. 


P = 0.925545352(94)). For the hyper-parameter length scale, 
we used a uniform prior based on the range of the instrumental 
parameter variations. For each instrumental parameter, we com- 
puted the maximum range of variations and set this as the max- 
imum possible value of the prior and the minimum was set to 
be one-sixth of this value to avoid over-fitting. For example, in 
the case of the roll angle, the prior is uniform between 60°and 
360°. We also included a jitter parameter for each light curve. 
The derived values for the jitter indicate that the errors provided 
by the DRP pipeline are slightly underestimated by a factor of 
approximately 1.2. 

For each transit, we started by using a 2D GP with the roll 
angle as the detrending instrumental parameter, then we anal- 
ysed the correlations of the residuals with the other instrumental 
parameters and added the instrumental parameter with the high- 
est correlations to a 3D GP. We performed model comparison to 
decide if the instrumental parameter should be added or not. In- 
strumental parameters were added if the difference of BIC was 
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Table 2. Priors for the fitted transit parameters. 


Parameter Prior Derived value 
To — Tree (days) U(-0.1;0.1) 0.000213 + 0.000062 
Rp/Rx N (0.1150, 0.0020) 0.11100 0.00014 
a/R, N(3.010, 0.013) 2.9829 + 0.0054 
inc [°] N(88.8, 1.1) 89.22 + 0.60 
LD1 N(0.5269, 0.0218) 0.5269 + 0.0218 
LD2 N(0.1279, 0.046) 0.1279 + 0.046 
Notes. 


T,,; = 2456836.2963007, U(a; b) is a uniform distribution between a 
and b; N (a; b) is a normal distribution with mean a and standard devia- 
tion b. 


higher than 3, indicating positive evidence in favour of the more 
complex model. The procedure was repeated until the residuals 
showed correlations with the instrumental parameters less than 
596 or they were already tested with model comparison. The first 
three light curves were highly correlated with the background 
and hence we tested a 2D GP with only the background as a 
detrending instrumental parameter. This was found to be better 
only for light curve number 3. This procedure was applied to 
both the OPTIMAL and the DEFAULT apertures and we chose 
the one with smaller final residuals. The instrumental parameters 
chosen with this procedure for the decorrelation of each of the 
transit light curves as well as the aperture chosen are given in 
Table [I] 


For parameter inference, we used the affine-invariant 
Markov-chain Monte-Carlo ensemble sampler implemented in 
emcee Foreman-Mackey et al.[2013). 
The fitting procedure was performed in two steps. First we per- 
formed a global fit using previously normalised light curves (first 
order polynomial based on out-of-transit data) and assuming a 
linear ephemerides and the best detrending GP model. From the 
global fit, we derived the transit parameters and their uncertain- 


ties similarly to[Barros et al.| (2020). These are given in Table [2] 


together with the priors used. In the second step, the posterior 
of the shape parameters a/R,, rp/R, and inc derived from the 
first step were used as priors for a second individual fit to each 
light curve to derive accurate transit times. In this second step, 
we simultaneously accounted for a linear normalisation of the 
transit parameterised by an out-of-transit level (Fou) and flux 
gradient (Fraa). Fitting the transit normalisation is important to 
avoid biases in the derived mid-transit times as shown in|Barros 
(2013). For each light curve, we used the best model GP 
determined in the previous step with the same hyper-parameters' 
priors mentioned above. The two step approach was adopted be- 
cause the correlations are different for each light curve and the 
detrending instrumental parameters considered are also differ- 
ent. Therefore, a simultaneous fit of the detrending instrumental 
parameters would imply a prohibitive number of parameters to 
fit (317). We performed tests that showed that this approach does 
not affect the results. 


The WASP-103b transit light curves together with the best 
multi-dimensional GP and transit model, chosen by model com- 
parison, are shown in Figure[3] The light curves show instrumen- 
tal effects that are well corrected by the GP model, as can be seen 
from the well behaved residuals. 


2.3.2. Pycheops 


For comparison, we also analysed the OPTIMAL extracted light 
curves with pycheops (Maxted et al.|2021), a custom python 
package developed specially for CHEOPS data. First we used the 
single visit analysis to determine the best parameters to use as 
decorrelation instrumental parameters. pycheops performs lin- 
ear decorrelation with several instrumental parameters such as 
contamination, background, position of the target on the CCD, 
and trigonometric functions of the roll angle and its harmonics. 
As in the multi-GP method, we used priors based on the transit 
parameters derived in[Delrez et al.] (2018). 

After analysing the data with the single visit model, we used 
the multivisit mode of pycheops to simultaneously fit the 12 
transits and determine the individual transit times. To avoid the 
large number of fitted parameters, pycheops has implemented 
a technique to perform an implicit decor- 
relation of several light curves using a GP. A detailed descrip- 
tion of pycheops and example applications to CHEOPS data are 
given in[Maxted et al.| (2021). The derived transit times with py- 
cheops are closer than | o to the ones derived with the multi- 
dimensional GP. The derived uncertainties in the transit times 
are also similar. 


3. Complementary observations 
3.1. Constraining the companion of WASP-103 


Using the lucky imaging technique, a possible companion to 
WASP-103 was detected at 0.242 + 0.016 arcsec by [Wóllert &| 
on 07 March 2015. They measured the posi- 
tion angle to be 132.7 + 2.7 degrees and contrast magnitudes 
to be Ai = 3.11 + 0.46 and Az = 2.59 + 0.35. These observa- 
tions were made with the AstraLux instrument (Hormuth et al. 
[2008). This companion was later confirmed by adaptive optics 
(AO) observations using the NIRC2 instrument at Keck 
on the 25 January 2016. They clearly detect a com- 
panion at 0.239 x 0.002 arcsec. Within the errors, no change of 
position was detected between the two observations. These ob- 
servations were used by |Cartier et al.|(2017) to perform a spectral 
energy distribution (SED) fit to the companion and estimate its 
parameters. This assumes that the companion is located at the 
same distance as WASP-103 which has not been confirmed. The 
stellar parameters derived by (2017) together with 
the position measurements derived by are re- 
produced in Table B] According to the AO measurements, if the 
companion is at the same distance as WASP-103A (552 +33 par- 
sec — Gaia EDR3 (Gaia Collaboration et al.[2021) - Table j), it 
would be at 131.9 + 8au from WASP-103. If the orbit is circu- 
lar, this would imply a period of 1114 + 103 years and a radial 
velocity (RV) signature with an amplitude of ~ 1334 m/s. 

There is an excess of astrometric noise in the Gaia data that 
is consistent with the existence of a companion for this system. 
This noise was present in the DR2 data release and in the recent 
EDR3. Furthermore, the Gaia derived parallax changed from 
1.14 + 0.17 in DR2 to 1.82 + 0.11 in EDR3 which is a 3.40 
change, which could be due to a deviation from single-source 
behaviour induced by a companion. Therefore, the companion 
still seems to be close to WASP-103 at the present epoch. 


3.1.1. Lucky imaging observations 


To better constrain the companion of WASP-103 and confirm 
that it is bound, we performed new lucky imaging observations 
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Fig. 3. Left panel: Transit light curves of WASP-103b obtained with CHEOPS. We overplot our best fit model that includes a transit model and 
the GP model to account for systematics dependent on the instrumental parameters. Right panel: Residuals of the fit of the transit model and the 
GP model. For clarity the errors are only shown in the residuals. The light curves and residuals are offset vertically for clarity. 


Table 3. Parameters for the companion of WASP-103A. 


Parameter Value and uncertainty 


Effective temperature Ter [K] 4330 + 100 
Surface gravity log g [g cm^?] 4.604 + 0.016 
Spectral type K5V 
Stellar mass M, [Mo] 0.721 + 0.024 
Stellar radius relative to A Rg/Ra 0.52 +0.05 
Distance to WASP-103 A [mas] 240.0 + 1.5 
Distance to WASP-103 A [au] 131.9 8 
Position angle PA [degrees] 131.37 0.35 


Notes. 

The stellar parameters of the companion of WASP-103 were derived 
from a SED fit assuming the same distance as WASP-103A 
[et_al.[2017). The position of the companion relative to WASP-103A 


was derived by (2016) . 


of WASP-103. We used the same instrument that discovered the 


companion (Wóllert & Brandner|2015), the AstraLux camera 
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(Hormuth et al.|2008) mounted on the 2.2-metre telescope at 


Calar Alto observatory in Almería, Spain. The observations were 
performed under excellent conditions (seeing 0.7 arcsec) in two 
filters SDSSi and SDSSz. We obtained 90,000 frames, each with 
an exposure time of 10 ms. 


As shown in Figure[A.1]and[A 2]of the appendix, we did not 
recover the companion found by|Wóllert & Brandner|(2015) and 
subsequently characterised 


(2017). We computed the sensitivity limits for our images by 
using the approach described in [Lillo-Box et al.| 2012] 2014). 
According to our contrast curve, we should have detected the 
companion if it was in the same location (separations and po- 
sition angle). No difference in the position of the companion 
was detected between the two previous observations 
[Ngo et al.|2016), but they were only sepa- 
rated by 10 months. Assuming that the target observed by us 
and the previous publications is the same, in order to explain 
the non-detection of the companion in our AstraLux images, the 
companion should have moved towards WASP-103A by at least 
~ 0.14 arcsecs, which corresponds to ~ 77 au at the distance to 
WASP-103 (TableH). 
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This is difficult to explain in a scenario where the companion 
is bound to WASP-103 since the time difference between the 
observations is just 6 years. The maximum velocity of a bound 
object is lower than the velocity at the periapsis given by the 
following: 


v2 = 2G(M. + Mceomp) 


i 2 
max d ( ) 
(Murray & Dermott|1999), where G is the gravitational constant 


and d is the distance of the periapsis. Assuming the periapsis 
is equal to 131.9au (Table B). in 6 years, the upper limit on the 
distance travelled by a bound object is 6.4 + 0.2 au. Therefore, 
we conclude that either the star is not bound (and hence we are 
seeing the relative proper motion of both stars, with the back- 
ground star disappearing behind WASP-103) or much less likely, 
unknown systematics have prevented its detection in the new ob- 
servations. Further high resolution images of WASP-103 or the 
Gaia DR3 will shed light on this system. 


3.1.2. CORALIE RV observations 


If WASP-103A has a bound stellar companion, we expect a long- 
term variation in the observed RVs. Previous RV observations 
of WASP-103 from the discovery paper and a follow-up paper 
do not show any long-term 
variation, but they only span 450 days. To constrain the longer- 
period RV variations, we obtained new RVs of WASP-103 with 
CORALIE, the same instrument that was used in previous stud- 
ies. Ten new observations were made between 18 March 2021 
and 30 May 2021 which increased the time span of the observa- 
tions to 8 years. 

The CORALIE spectrograph is installed at the Nasmyth fo- 
cus of the Swiss Euler 1.2m telescope and 
has a spectral resolution of 60,000. The light can be injected 
through two fibres allowing it to observe the science target and 
to perform a simultaneous monitoring of the sky or a wavelength 
calibration with a Fabry Perot etalon. In November 2014, the 
spectrograph benefited from a major upgrade, which introduced 
an RV offset which we modelled as a simple offset between the 
two sets of data. The observations were processed with the stan- 
dard data reduction pipeline. The RVs were derived with the 
weighted cross-correlation technique and a 
G2 mask was used as it is optimised for late F-type stars such as 
WASP-103. 

The RVs were analysed with the code LISA 
2021) which uses the radvel python package (Ful- 
ton et al./2018) to model the RV observations. The RV model is 
parameterised by the semi-amplitude of the RV signal (K), the 
planetary period (P), the mid-transit time (To), and the products 
of the planetary eccentricity by the cosine and sine of the stellar 
argument of periastron e cos w, esin c. We used Gaussian pri- 
ors for the planetary orbital period and mid-transit time based 
on the constant period model derived in Section[4] We included 
an offset between the new RV observations and the previously 
published RV observations to account for the RV shift due to the 
upgrade of the instrument mentioned above. We also included 
a jitter parameter for each dataset to account for unknown sys- 
tematic noise or short-term stellar activity. We compared a model 
with a beta prior on the orbital eccentricity with 
a circular model. Due to the short time span of the observations 
compared with the expected orbital period of the possible binary, 
the visual companion, if bound, would lead to a trend in the RV 
observations. Moreover, given the short span of the two epochs 
of observations and the large gap between them, this trend can be 


represented by an offset between both observations. Therefore, 
the fitted offset is a combination of the instrumental offset and 
the trend due to the possible companion star. 

With the free eccentricity model, we found a non-significant 
eccentricity of 0.11 + 0.06. The difference between the Bayesian 
information criterion (BIC) of the eccentric and circular model 
is 0.011 which implies that the eccentric model is not justi- 
fied. As mentioned above, due to the very short orbital period 
of WASP-103b, the orbit of the planet is expected to be circu- 
larised and the rotation of the planet synchronised with the or- 
bital period. Therefore, we adopted the circular model for the 
planet. We found a semi-amplitude K = 268 + 14 m/s in agree- 
ment with previous results (Delrez et al./2018). We also found 
an offset between the previous observations and the new obser- 
vations of 14 + 45 m/s. At 3c we can exclude an offset higher 
than 151 m/s and lower than —119 m/s. The relative offset due 
to the instrumental upgrade between the two observations is ex- 
pected to be between 14 and 24 m/s (CORALIE team private 
communication). Therefore, we conclude that there is no signif- 
icant offset between the new and the previous observations of 
WASP-103b. At 3a we can also exclude RV variations higher 
than 151 — 24 = 127 m/s and lower than —119 — 24 = —143 m/s 
over 8 years. As outlined in Section [4.1] this limit on the ampli- 
tude of the RV variation does not allow us to discard the bound 
scenario. 


3.2. Stellar parameters 


Thanks to the new data release of Gaia, the stellar parame- 
ters of WASP-103A can be better constrained, which in turn 
allowed us to better constrain the mass and radius of the exo- 
planet WASP-103b. Using a modified version of the infrared flux 
method (IRFM; [Blackwell & Shallis/1977), we determined the 


radius of WASP-103A via relationships between various stellar 


parameters recently detailed in |Schanche et al.|(2020). We con- 


structed the SED from stellar atmospheric models using the stel- 
lar parameters from SWEET-Cat as priors, 
and subsequently attenuated the SED to account for reddening. 
The SED was further corrected for the companion using the cal- 
culated contamination estimate from the stellar parameters for 
the companion in |Cartier et al.|(2017) and reproduced in Table[3] 
The corrected SED was then convolved with broadband response 
functions for the chosen bandpasses to obtain synthetic photom- 
etry which allowed us to compute the bolometric flux, and hence 
the radius, of the target. We retrieved broadband fluxes and un- 
certainties from the most recent data releases for the following 
bandpasses: Gaia G, Ggp, and Grp; 2MASS J, H, and K; and 


WISE W1 and W2 (Skrutskie et al.|2006| Wright et al.|2010 
Gaia Collaboration et al.|2021). We also used the arras cata- 
logues (Castelli & Kurucz|2003) of model stellar SEDs. Within 


the IRFM, the distance used to convert the angular diameter of 
WASP-103A to the stellar radius was calculated from the Gaia 
EDR3 parallax with the parallax offset of|Lindegren et aT. (2021) 
being applied. Using a Markov chain Monte Carlo (MCMC) fit- 
ting approach, we estimated the stellar radius of WASP-103A to 
be R, = 1.716 + 0.119 Re. This is larger than the previous esti- 
mates, that is 1.436 + 0.052 (Gillon et al.|2014) due to the greater 
distance to WASP-103 derived from the EDR3 parallax. 

We derived both the stellar mass M, and age f, by em- 
ploying two different sets of stellar evolutionary models, namely 


PARSEC v1.2! |(Marigo et al.[2017) and CLES (Code Liégeois 


! PAdova & TRieste Stellar Evolutionary Code: http://stev. 
oapd.inaf.it/cgi-bin/cmd 
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d'Évolution Stellaire 2008). The input parame- 
[ 


ters we used were the stellar [Fe/H], Ter, and Rẹ to locate the 
star on the Hertzsprung-Russel Diagram (HRD) plus the pro- 
jected rotational velocity v sin i, which was plugged into the gy- 
rochronological relation by [Barnes] (2010) to remove isochronal 
degeneracies; for further details, see 
82.2.3). We performed a direct interpolation of the input param- 
eters within pre-computed grids of PARSEC models thanks to 
the isochrone placement algorithm presented in [Bonfanti et al.| 
and {Bonfanti et al.| (2016), obtaining the first pair of age 
and mass values. The second pair was inferred by directly com- 
puting the evolutionary track through the CLES code and then 
choosing the best-fit solution following the Levenberg-Marquadt 


minimisation scheme (Salmon et al.|2021). After checking the 


consistency of the two mass and age values following the y? 


criterion discussed in detail in |Bonfanti et al.| (2021), we com- 


bined the respective distributions of age and mass together to 
obtain our final estimates of M, and t,. The mass is in agree- 
ment with previous estimates, while the age is much better con- 


strained (Gillon et al./2014 2018). The final stellar 


parameters and the 1 c uncertainties are given in Table [4] 


Table 4. Stellar parameters of WASP-103A. 


Parameter Value and uncertainty 


Effective temperature Teg [K] 6013 + 44 
Surface gravity log g [g cm?] 4.24 + 0.15 
Iron abundance [Fe/H] [dex] 0.08 + 0.04 
Spectral type F8V 
Parallax* p [mas] 1.8110 + 0.1073 
Distance to Earth d [pc] 552 + 33 
Stellar mass M, [Mo] 1.204 + 0.046 
Stellar radius R, [Ro] 1.716 + 0.119 
Stellar age t [Gyr] 5.2 + 0.8 
Stellar luminosity L, [Lo] 3.47+ 0.49 


Notes. 


*Parallax from Gaia EDR3 (Gaia Collaboration et al.|2021) using the 
formulation of (2021). 


3.3. Re-analysis of previous transits 


Given the very high signal-to-noise required to detect the tidal 
deformation, we have re-analysed high signal-to-noise transits 
of WASP-103 previously obtained with the Spitzer and Hubble 
space telescopes. These are combined with the 12 new transits 
obtained with CHEOPS in our final analysis. 


3.3.1. Spitzer observations 


We re-analysed the Spitzer archival data of WASP-103b which 


has already been published (Kreidberg et al.|2018). We down- 
loaded WASP-103b archival IRAC data from the Spitzer Her- 


itage Archive (http://sha.ipac.caltech.edu). The data 
consist of one full phase curve of WASP-103b at 4.5 um (chan- 
nel 2) and one at 3.6 um (channel 1), both were obtained under 
program ID 11099 (PI L. Kreidberg) taken on 19 and 28 May 
2015. The reduction and analysis of these datasets are similar to 
(2016a). We modelled the IRAC intra-pixel sensi- 
tivity (Ingalls et al.|2016) using a modified implementation of the 
BiLinearly-Interpolated Sub-pixel Sensitivity (BLISS) mapping 


algorithm (Stevenson et al.|2012). We used a modified version 
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of the BLISS mapping (BM) approach to mitigate the correlated 
noise associated with intra-pixel sensitivity. In our photometric 
baseline model, we complement the BM correction with a lin- 
ear function of the point response function (PRF) full width at 
half-maximum (FWHM). 

In addition to the BLISS mapping, our baseline model in- 
cludes the PRF’s FWHM along the x and y axes, which signif- 
icantly reduces the level of correlated noise as shown in previ- 


ous studies (e.g. 2014 20 16a; |De-| 
; (Gillon et al. Mendonga et al. 


Our baseline model does not include time-dependent parame- 
ters. Our implementation of this baseline model is included in 
an MCMC framework already presented in the literature 
et al.12012). We ran two chains of 200,000 steps each to deter- 
mine the phase-curve properties and obtained the best detrended 
transit light curves which are analysed together with HST and 
CHEOPS transits in Section] From our BM+FWHM baseline 
model, we obtained a median root mean square (RMS) of 3450 
ppm per 10.4s integration time at 3.6 wm and 4480 ppm with the 
same integration time at 4.5m. 


3.3.2. HST observations 


We re-reduced the Hubble transit observations taken on 26-27 
February 2015 and 2 August 2015 with HST Program 14050 
which were originally published by [Kreidberg et al.|(2018). The 
target was acquired in both forward and backward scanning di- 
rection using an exposure time of 103 s. We used the frames in 
the IMA format, each one containing 16 non-destructive reads 
(NDR; Messer which were pre-processed by the 
CALWEC3 pipelind*| version 3.5.2. 

Wavelength calibration, NDR operations, background sub- 
traction, cosmic ray and bad pixels rejection, and correction for 
drifts were carried out following standard procedures, as de- 
scribed in [Bruno et al.|(2018] and references therein). Then, we 
integrated the stellar spectra in the 1.115-1.625 um wavelength 
range to obtain the band-integrated transits. Following standard 
practice (e.g. Deming et al.[2013), we rejected the first HST or- 
bit of the transit obtained on 2 August 2015, which was at the 
beginning of the phase curve observation, and the first data point 
of every orbit for both transits. 


We then used a method similar to|Kreidberg et al.| (2018) to 


remove the instrument systematics with the model described in 


(2014), that is with a second-order polynomial 


and an exponential ramp, 


S(t) = CA +1790  ri))(1— e?9*^ + rao), (3) 
where we fitted for C and ro.4, and 0 and ¢ represent the plane- 
tary and HST phase, respectively. It was also necessary to add a 
shift to the HST orbital phase, ø = 2z[(t—V) mod Pys7]/Pusr, 
where y = —0.045 d is for the February 2015 visit and y = 
—0.1 d is for the August 2015 visit, respectively, and Pysr is 
Hubble's orbital period. 

The systematics were fitted simultaneously with the transit 
model of a non-deformed planet using the model of 


& Agol| (2002) (implemented in |Kreidberg| (2015b) software) 
and scipy's optimize.minimize function (Virtanen et al.|2020 


and references therein). The best detrended transit light curves 
are analysed together with Spitzer and CHEOPS transits in Sec- 
tion [5|The mid-times of each exposure were converted to BJD- 


? https://www.stsci.edu/hst/instrumentation/wfc3/software- 
tools/pipeline 


S. C. C. Barros et al.: Detection of the tidal deformation of WASP-103b at 3 o with CHEOPS 


TDB using the online applet based on the method of 


fal 2010). 


4. Period evolution of WASP-103b 


Using the procedure outlined in Section [2.3.1] we obtained the 
mid-transit times of the CHEOPS observations. These are given 
in Table [5] In this table we also included the mid-transit times 
derived for the Spitzer and HST transits. We reiterate that the 
derived CHEOPS transit times obtained with our method are 
within 1 c of the ones derived using pycheops showing that our 
detrending methods are robust. 

We combined our derived mid-transit times (124-4) with the 
32 previously published mid-transit times of WASP-103 which 


were presented in Table 3 of [Maciejewski et al.| (2018), some 


of which are reanalyses of previously published values (Gillon 
et al.|2014] [Southworth et al.[2015} [Delrez et al.|2018] {Turner} 
et al.|2017] |Lendl et al.|2017). We also added the four transit 
times Mise aL, in Patra et al.| (2020 (2020). Therefore, 
in total we have 52 mid-transit times of WAGE 103b spanning 
seven years. 

For parameter inference, we used emcee as in Section 
We compared a linear ephemerides (constant period) model with 
a quadratic ephemerides (constant derivative period) model. We 
included a multiplicative jitter parameter in our analysis. 

For the constant period model, we obtained Tọ = 
2457511.944458*0 00" (BJDrpg) and P = 0.925545485 + 
0.000000049 days. The BIC of this fit is 79.8. We found a jit- 
ter of 1.18. 

We considered a constant a ae model with the 


following form (e.g. Maciejewski et al 020): 


. E(E-1) 
Tui = To + PX E+ PPX a 


(4) 
where E is the transit epoch and P is the period derivative. 

For this model we derived Tg = 2457511.944344 + 
0.000075(BJD rpg), P = 0.9255453 + 0.000000089 days, P = 
3.5 + 1.8 x 1071? days/day, and jitter = 1.15. The jitter is slightly 
lower than for the linear model. We found a hint of an increas- 
ing orbital period, contrary to what was expected if tidal decay 
was dominating the orbital evolution of the system. The period 
derivative was found to be positive at 2.1 c which is not signif- 
icant. The BIC of the quadratic model is 78.05, giving a differ- 
ence of BIC in favour of the quadratic model of only 1.8. There- 
fore, according to the BIC, the added complexity of the quadratic 
model is not strongly justified and the linear ephemerides is pre- 
ferred. 

Under assumption that the period variation observed is due 
to tidal decay (i.e. the period is actually decreasing and the varia- 
tion seen is due to statistical uncertainties), we can derive a lower 
limit to the tidal dissipation parameter using the following equa- 


tion (e.g. [Maciejewski et al.|2020): 


/ 


*' 2 M. 


272 M, (RV 1 
an *( ) 5 (5) 


where M, and M. are the planetary and stellar masses, respec- 
tively, R. is the stellar radius, and a is the semi-major axis of the 
planet's orbit. We derived a lower limit on the period derivative 
to be —1.3 x 107!° days/day at the 99.7% confidence interval. 
This implies that the tidal dissipation parameter is higher than 
1.6 x 10° at 3c, corresponding to a 99.7% confidence interval 
if we assume that only tidal decay affects the period derivative. 


IR 
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Fig. 4. Derived mid-transit times of WASP-103b after removing a lin- 
ear ephemerides. The CHEOPS data are shown in blue, the re-reduced 


HST data and Spitzer are shown in magenta, and the previously pub- 
lished times are s D black. The previously published quadratic 


ephemerides (Patra et al.12020) is shown in green with the 1 c uncer- 


tainty limit, while our new esti is shown in red. 


This limit is more than an order of magnitude higher than pre- 
vious studies that found Q^ > 1.1 x 10? at 95% 


for WASP-103b. At 9596 confidence, our results allow us 
to exclude a negative period derivative. The quadratic fit to the 
derived transit times is given in Figure |4| We found a period 
derivative that is smaller than the previous estimation 
[2020). although the higher precision results in a higher signifi- 
cance for being positive. 

Figure 4| shows that the first two CHEOPS transits have a 
slightly late mid-transit time compared with the other observed 
CHEOPS transit times, although consistent within the errors. 
This is probably due to the difficulty in detrending CHEOPS data 
when the duration of the visit is shorter than three CHEOPS or- 
bits. It is also known that transits with poorly covered ingress or 
egress can lead to biases in the derived transit times 
2013). To test if this could influence our results we repeated the 
linear and quadratic model fits excluding the first two transits. 
We found no significant differences in the derived model param- 
eters or model comparison. We also repeated the fits using the 
transit times derived with pycheops instead of the ones derived 
with the multi-dimensional GP and found the same results. 

Although it is likely that the observed period variation is due 
to statistical dispersion and that the orbit is decaying due to tides, 
it is interesting to explore other factors that could affect the or- 
bital evolution of this system. In the next subsections, we explore 
scenarios that could explain an increase in the orbital period in 
case it becomes significant after future observations. 


4.1. RV acceleration due to a companion 


The existence of a companion of WASP-103 could lead to 
RV acceleration which would produce transit timing variations 
due to a change in the light travel time. Assuming a quadratic 
ephemerides and that the observed period derivative is due to 
the Doppler effect of a line-of-sight acceleration (Pay), we can 
derive the line-of-sight acceleration (a,) using the following: 


, a,P 
Pry = , 
c 


(6) 


where c is the speed of light. We obtained a, = 0.113 + 0.058 
m/s/day. 
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Table 5. Derived mid-transit times of WASP-103b. 


Epoch time T4 g- 
BJD TDB -2450000 (days) (days) (days) 
-466 77080.64040 0.00033 0.00028 
-377 7163.01442 0.00036 0.00036 
-368 7171.34309 0.00024 0.00024 
-297 7237.05769 0.00027 0.00030 
1563 8958.57247 0.00037 0.00035 
1578 8972.45609 0.00024 0.00025 
1581 8975.23171 0.00036 0.00034 
1593 8986.33822 0.00039 0.00031 
1596 8989.11503 0.00024 0.00020 
1603 8995.59388 0.00022 0.00023 
1616 9007.62623 0.00023 0.00022 
1617 9008.55119 0.00022 0.00023 
1624 9015.03013 0.00022 0.00022 
1629 9019.65833 0.00023 0.00023 
1630 9020.58337 0.00025 0.00023 
1631 9021.50904 0.00025 0.00029 
Notes. 


The first four entries refer to the HST and Spitzer transits, while the last 
12 were derived from the new CHEOPS observations. 


During the eight years since the discovery of WASP-103b 
until now, this would imply an RV variation of 330 + 168 m/s. 
This is higher than the RV offset that we measured in Sec- 
tion[3.1.2] but it is still compatible within the errors. 

We can also calculate the expected acceleration from the vi- 
sual companion of WASP-103 if it is bound (discussed in Sec- 
tion|3.1) using Newton's law: 


GM, comp 
a= ; 
d? 


(7) 


where Mcomp is the mass of the visual companion star and d is 
the separation between the two stars. In this case, where the com- 
panion is observed to be spatially separated from WASP-103, we 
need to account for the projection of the acceleration in the line- 
of-sight given that only the radial component of the acceleration 
results in a variation of the observed period. If the angle between 
the line-of-sight and the companion is 6, then a,44 = a cos 0 and 
d = ó[ sin 0. Where 6 is the projected separation of the star and 
the companion that we previously derived (SectionB.1] as 134+8 
au. Hence, 


comp 


c o £2 
ärad = —g Cin 0)* cos 0. 


(8) 


The maximum value of equation [8] is obtained for cos@ = 
P Assuming Mcomp = 0.721 + 0.024 Mo(Table 3). we can set 
an upper limit on a,44 < 0.00796 + 0.00095 m/s/day. This cor- 
responds to an RV acceleration of 23.6 + 2.7 m/s in 8 years. 
This is compatible with the RV offset that we measured between 
the new CORALIE observations and the previously published 
observations (Section (3.1.2). The difference between the accel- 
eration expected if the transit timing variations are due to ac- 
celeration from a bound companion and the acceleration from 
the observed visual companion is 0.105 + 0.057 m/s/day. Hence, 
the acceleration produced by the visual companion is probably 
not enough to produce the period derivative estimated with the 
quadratic model. However, we cannot exclude it at more than 
2:01 
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Long-term RV monitoring of WASP-103b and the next Gaia 
DR3 will allow us to better constrain the existence of possi- 
ble bound companions to WASP-103 and correct the line-of- 
sight acceleration light travel time, allowing us to better con- 
strain the tidal dissipation parameter. The hypothesis that the vi- 
sual companion observed by lucky imaging and AO observation 
is responsible for the transit timing variations and the offset in 
CORALIE RVs cannot be completely rejected. However, the ab- 
sence of this star in our new lucky imaging observation and the 
fact that the predicted acceleration by this star is 2c lower than 
required to match the observations, suggests that other mech- 
anisms would be required to explain an increase in the orbital 
period of the planet. 


4.2. Applegate effect 


Eclipse times of binary stars have been shown to vary due to 
variations in the quadrupole moment of the stars driven by stellar 
activity. Low mass stars with convective outer layers have varia- 
tions of their quadrupole moment due to a distribution of angular 
momentum driven by stellar activity cycles. The change of the 
quadrupole moment of the star leads to quasi-periodic variations 
of the eclipse times of the companion over timescales of years 
to decades. This effect has been measured in many eclipsing bi- 


naries and is known as the Applegate effect (Applegate| 1992). 


Another explanation for the observed period changes in binaries 
was proposed by |Lanza et al.|(1998). In this case, the Applegate 
effect would be due to a cyclic transformation of rotational ki- 
netic energy into magnetic energy and back to rotational kinetic 
energy. If the Applegate effect is detected, it would allow us to 
probe the nature of the dynamo mechanism of low mass stars 


(Lanza et al.[1998). 

Since exoplanets host low mass stars with some dynamo ac- 
tivity, it is expected that the Applegate effect is also present in ex- 
oplanet systems; although, it has never been observed. [Watson &| 
estimated the transit time variations due to the Ap- 
plegate effect for a few transiting exoplanets. They show that for 
stellar dynamos with timescales of 11 years, the Applegate effect 
is less than 5 seconds for most exoplanet host stars. However, 
for stellar dynamos with longer timescales, the effect can reach 
a few minutes. Using their equation 13 and assuming the stel- 
lar parameters given in Table [4] the semi-major axis of the orbit 
a = 0.01985 au, the observation time span T = 7 years, and es- 
timating the stellar angular rotation velocity of WASP-103 from 
the vsini given in Gillon et al] 014) (10.6 kms~!), we con- 
clude that the Applegate effect in WASP-103b would produce 
transit timing variations < 38 seconds over the time span of the 
available observations. Assuming the quadratic ephemerides, we 
found that at the mid-epoch of CHEOPS observations, the mea- 
sured transit time, is 1.69 + 0.81 minutes later than what would 
be expected by a linear ephemerides. Therefore, this is higher 
than the expected timing variations from the Applegate effect, 
although in agreement at 1.3 o. Hence we conclude that the Ap- 
plegate effect could be affecting the measured transit times of 
WASP-103b. 


4.3. Apsidal precession 


If a planet’s orbit is slightly eccentric, then its orbit would be 
apsidally precessing. For hot Jupiters, the precession timescale 
is expected to be decades. In this case, there is a long-term os- 
cillation of the apparent period. Modelling the period variation 
would allow us to determine the planet Love number and con- 
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strain its internal structure. For WASP-103b, a zero eccentricity 
was assumed by and favoured by the analysis 
of (2018) and our own analysis using the new RV 
measurements presented in Section [3.1.2] 

Nevertheless we attempted to fit the times of WASP-103 with 
a transit timing model assuming apsidal precession (Patra et al. 


ROT: 


eP, 
T mia = To + Ps X E- — x cos w, (9) 
nz 
where, 
dw/dE 
Be a(1- p l (10) 
2n 


We found that the apsidal precession model is a poor fit to the 
transit times with a BIC = 82.7. This is due to the two extra 
parameters compared with the quadratic model that is already 
not justified by the BIC compared with the linear model. How- 
ever, we obtained physical values for the fitted parameters. We 
obtained 72 = 1.10412, x 107? rad, e = 0.00054*90055 <, and 
w = 0.49107 rad. 

The most important contributions to the apsidal precession 
rate of hot Jupiters are those coming from the tides raised by the 


star on the planet and from the rotation of the planet 
2009). Assuming synchronous rotation, the leading term 


in the expression of this precession rate at low eccentricity is as 


follows (see, e.g. eqs (6)+(10) of|Ragozzine & Wolf|2009): 


5 
E31 = l6z(hy — Du (=) , 

T+R p\4a 
Using our fitted value for the precession rate of 1.10412, x 107° 
rad, we can estimate the Love number to be hy = 1.35 + 0.43 
which is compatible with our estimate from the deformation of 
the light curve (Section [5}. 

Current observational constraints on the eccentricity cannot 
rule out such a small value ~ 0.00054. However, due to the 
short circularisation timescale, the eccentricity of WASP-103b 
is expected to be zero unless there is an external perturbation. 
For example, a planetary companion can excite the eccentricity 
of WASP-103b. The eccentricity can also be excited by gravi- 
tational perturbations from the star’s convective eddies as pro- 
posed by Phinney (1992). 

Therefore, we cannot completely rule out that apsidal pre- 
cession is affecting the transit times of WASP-103b given our 
current constraints on the eccentricity of the planet although this 
is, a priori, not expected. Future monitoring of the transit times of 
WASP-103b can disentangle apsidal precession from tidal decay 
since for apsidal precession the variations are sinusoidal. The 
times of occultation can also be used to disentangle both sce- 
narios because in the apsidal precesion, these are anti-correlated 
with the times of transit. 


(1) 


5. Tidal deformation analysis 


As mentioned above, WASP-103b is the exoplanet with the high- 
est expected deformation signature due to its large radius and 
close proximity to its host star. We attempted to measure the 
deformation and tidal Love number of WASP-103b, combin- 
ing the 12 new high-precision transits obtained with CHEOPS 
with two HST transits and two Spitzer transits (3.4 and 4.5 mi- 
crons). To model the tidal deformation, we used the implemen- 


tation of Akinsanmi et al.| (2019) based on the parametrisation 


Table 6. Priors for the fitted transit parameters. 


Parameter Prior 

a/R, U(2.5, 3.5) 

b U0, 1) 

hy U(0, 2.5) 

log On N(-6.7581, 0.0534) * 
R,/R, each instrument U(0.05, 1.5) 


c CHEOPS N(0.7045, 0.0147) 
a CHEOPS N(0.7670, 0.0199) 
c HST N(0.5714, 0.0218) 
a HST N(0.4285, 0.0175) 
c Spitzer 1 N(0.2772, 0.0085) 
æ Spitzer 1 (0.4730, 0.0229) 
c Spitzer 2 N(0.2280, 0.0067) 
a Spitzer 2 N (0.483, 0.021) 
Notes. 


* Derived from the RV analysis in Section 


of (2014). This implementation uses the ellc transit tool 
(Maxted|2016) and it is also freely available. The model param- 
eters are the normalised separation of the planet (a/R,), the im- 


pact parameter (b), the Love number (hp), the logarithm of the 
planet-to-star mass ratio multiplied by the sine of the inclination 
(In Q,, = In (Iz sin inc)), and, for each filter, the planet-to-star 
radius ratio (Ry/R,) and the power-2 limb darkening (LD) coef- 
ficients (c and o). Following equation[I] in this ellipsoidal model, 
the radius of the planet is parameterised by the volumetric radius 
Ry = Wriror3. The LD coefficients were parameterised accord- 
ing to and [Short et al.] (2019) to minimise the 
correlations between them and to avoid non-physical solutions. 

The priors for each parameter are given in Table [6] For the 
shape parameters, we used uniform uninformative priors instead 
of normal distributions based on previous data because the pre- 
vious results were obtained assuming sphericity, which impacts 
the derived shape parameters. We assumed the period and mid- 
transit times derived in Section 4] We included a multiplicative 
jitter term for CHEOPS, HST, and Spitzer channel 1 and chan- 
nel 2 to account for any underestimation of the uncertainties. 
For each light curve, we corrected the contamination due to the 
visual companion star (see Section[3.1), assuming the stellar pa- 
rameters given in Table [3] and the respective filter transmission 
functions. 

For the parameter inference, we used the nested sampling al- 


gorithm implemented in Dynesty (Speagle|2020 
2019} |Skilling|2012, |2004) which provides posterior estimates 


and also the Bayesian evidence useful for model comparison. We 
fitted the tidal deformation parameterised by ^; and compared it 
with a spherical model (^; = 0). The comparison of the models 
illustrates biases in the derived shape parameters if the defor- 
mation is ignored. The derived transit parameters are given in 
Table [7] for the spherical and the ellipsoidal model. The derived 
jitter parameters for the ellipsoidal model are 1.00, 1.11, 1.21, 
and 1.09 for Spitzer channel 2, Spitzer channel 1, CHEOPS, and 
HST, respectively, showing that our errors are robust and the de- 
trending was successful. The derived jitter parameters for the 
spherical model are similar to the ellipsoidal model. 

We overplotted the best model that accounts for tidal defor- 
mation on the time-folded light curves of WASP-103b taken with 
Spitzer 2, Spitzer 1, CHEOPS, and HST in Figure [5] We also 
show the residuals of the spherical model and overplotted the 
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Table 7. Derived transit parameters of WASP-103b for an ellipsoidal planet model and a spherical model. 


Parameter Spherical Model fit (S)  Ellipsoidal Model fit (E) 
Love number hy - 1.59108 
In Qm - —6.761 + 0.049 
40.0061 40.0046 
s par aaa PRI" 
OY 0.039 DUET 0.027 
M,/Miup 1.464 + 0.096 1.460 + 0.089 


Ry/R, Spitzer 2 
Ry/R, Spitzer 1 
Ry/R, HST 

Ry/R, CHEOPS 


0.12245 + 0.00086 
0.11928 + 0.00059 
0.11639 + 0.00018 
0.11588 + 0.00020 


0.1297 + 0.0028 
0.1257 + 0.0024 
0.1222 + 0.0021 
0.1215 + 0.0020 


Ry [Rjup] Spitzer 2 2.044 + 0.142 2.165 + 0.155 
Ry [Ryup] Spitzer 1 1.992 + 0.138 2.097 + 0.150 
Ry [Rjup] HST 1.943 + 0.135 2.039 + 0.144 
Ry [Rjup] CHEOPS 1.935 + 0.134 2.027 + 0.145 
Pp [pup] Spitzer 1 0.171200 0.144200 
Pp [PJup] Spitzer 1 0.1850 935 0.15875 934 
Pp [sup] HST 0.199 Fg 0.173994 
Pp [P.sup]| CHEOPS 0.202*9050 0.175700 


difference between the best fit spherical model and the best fit 
ellipsoidal model. As shown by (2019), this 
is the measurable signature of the deformation of a planet in a 
transit light curve. This signature has two components. The first 
one is the signature of the oblateness (r2 > r3 ) resulting in an 
oscillation in the residuals of the flux during ingress and egress 
(Seager & Hui[2002] Barnes & Foriney 2003). The second one 
rises from rı > r2 due to the change of the projected area of the 
ellipsoidal planet as it rotates synchronously with its orbit. This 
results in a bump that has its maximum at the minimum of the 
projection which is the middle of the transit (Correia]2014). A 
schematic view of the geometry of how the deformation changes 
a transit light curve is given in Figure A.1 of|Correia|(2014). The 
change in the amplitude of signature of the deformation with the 
wavelength of the observations due to the change in the limb 
darkening and the larger planetary radius at longer wavelengths, 
as it can be seen in Figure |5| is noteworthy. This prevented us 
from phase folding all light curves and signatures so that our 
results could be visualised better. 


In Figure [6] we show the correlation plots and the poste- 
rior probability distributions for the derived transit parameters of 
WASP-103b. As expected, there is a large correlation between 
the Love number and the radius ratio for the ellipsoidal model 
that leads to a larger uncertainty of the parameters of this model. 
For simplicity, we do not show the distribution of the LD param- 
eters and the jitter parameters because their shape is very well 
approximated by a Gaussian and they are very similar for the 
two models. 


We derived the radial Love number of WASP-103b to be 
hy = pu. This is the first time that a 3 o detection of the 
Love number has been achieved for an exoplanet directly from 
the analysis of the deformation of the transit light curve. To ob- 
tain this result we combined the datasets from the three instru- 
ments: CHEOPS, HST, and Spitzer. To show the importance of 
each data set, we fitted the Spitzer and HST light curves sepa- 
rately and together. These results are given in Table[8]and show 
that the addition of the CHEOPS data was necessary to obtain a 
3c detection. It also justifies that the signature is not evident by 
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Table 8. Comparison of the derived Love number for the individual 
instruments and their combination. 


Data set Love number Significance Bayes Factor 

SP2, SPI 1.362579 1.70 1.7 

HST 0.99755 1.70 0.71 

CHEOPS 1.7475 49 2.50 6.1 

HST, SP2, SP1 L16 1.8 o 1.0 

All data 1.594553 3.00 9.1 (17*) 
Notes. 


* Corrected value of the Bayes factor as explained in Section|5.1] 


eye in Figure[5]for any individual datasets. We conclude that the 
Love number of WASP-103b is similar to the Love number mea- 
sured for Jupiter (1.565 + 0.006 — [Durante et al.[2020), suggest- 
ing a similar internal structure despite the much larger radius and 
much higher levels of irradiation for this exoplanet. The derived 
Love number of WASP-103b is higher than the one estimated 
for HAT-P-13b by [Batygin et al.|(2009). This new measurement 


of the Love number can be used to lift the degeneracy of inter- 


nal composition models and allows the 
derivation of the core mass of WASP-103b similarly to what was 
done for HAT-P-13b Buhler et al.[2016). 
We found that the volumetric radius derived with the ellip- 
soidal model is 5-646 bigger than the radius estimated with the 
spherical model. Therefore, not accounting for deformation bi- 
ases the derived planetary radius and hence the planetary density 
(^ 14%) and composition. This is the first time that this bias that 
was predicted by and (2014) has 
been directly measured. The large tidal deformation in ultra-hot 
Jupiters affects their phase curve observations and consequently 
their atmospheric characterisation. Previous phase curve mea- 


surements of WASP-103b (Delrez et al.[2018| |Lendl et al.|2017] 
Kreidberg et al.|2018) have corrected tidal deformation using 
theoretical estimations (e.g. 2011} 2011) 


that assume an interior structure for the planet. Our measure- 
ment of the Love number will allow an assumption-free correc- 
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Ellipsoidal Model fit to all the datasets 
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Fig. 5. Time folded transit light curves of WASP-103b obtained with Spitzer 1, Spitzer 2, CHEOPS, and HST. The CHEOPS transit light curve is a 
combination of 12 individual transits, while the HST light curve is a combination of two transits. The best fit ellipsoidal transit model is shown in 
blue. We also plotted the residuals of the best fit ellipsoidal model (blue) and the best fit spherical model (red) binned to 5 minutes. On the latter, we 
overplotted the signature of the deformation (green) which is the difference between the best fit spherical model and the best fit ellipsoidal models. 
For clarity, we replotted a zoom of the signature of the deformation in the bottom panel for each filter. We also show the mean uncertainties of the 
original data points and of the binned residuals. 


tion based on direct observations. This will allow a more accu- account for the deformation of WASP-103b could affect the in- 
rate estimation of the day-side and night-side temperatures from — terpretation of its transmission spectra (Lendl et al.|2017). 
phase curve observations. It is also possible that neglecting to 
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Fig. 6. Derived correlation plots and posterior probability distributions of the transit parameters of WASP-103b for the spherical (red) and ellip- 
soidal model (blue). The vertical lines show the median of the distributions and the shaded area shows the 6896 confidence intervals. We show 
the 1 o (dark blue and dark red) and 2c (light blue and light red) contours. We obtained a 3 ø detection of the Love number. The parameter 
distributions also clearly show that the ellipsoidal model is not as well constrained as the spherical model due to strong correlations between the 
Love number and the radius ratio. For the ellipsoidal model, the radius ratio refers to the volumetric radius. The superscripts spl, sp2, ch, and hst 


refer to the two Spitzer channels, CHEOPS, and HST, respectively. 


5.1. Assessing the significance of the detection 


One way to assess the significance of the detection is to perform 
model comparison - probability of one hypothesis versus an- 
other. Bayesian model comparison requires computing the odds 
ratio between two hypotheses (e.g. 2014). The odds 
ratio is the multiplication between the prior odds and the Bayes 
factor (ratio of the Bayesian evidences). The prior odds are the a 
priori probability of each model. Given the strong tidal forces 
WASP-103b is subjected to by its host star, theoretically, we 
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know that the planet has to be deformed. Hence, the prior proba- 
bility of the spherical model is zero which implies that the odds 
ratio in favour of the ellipsoidal model is infinity and renders the 
Bayes factor irrelevant. Nevertheless, we estimated the Bayes 
factor of the ellipsoidal compared with the spherical model us- 
ing the evidence computed with Dynesty. We found the Bayes 
factor (ratio of the Bayesian evidences) of the ellipsoidal com- 
pared with the spherical model to be 9.1, giving positive evi- 


dence for the ellipsoidal model (Kass & Raftery|1995). However, 
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the Bayes factor penalises more complex models which is incor- 
rect in our case since, as mentioned above, the planet is expected 
to be deformed and not accounting for deformation significantly 
biases the derived transit parameters, especially the planetary ra- 
dius. To correct the penalisation of the extra parameters, we fitted 
an ellipsoidal model with a fixed value of A; and In Qm, corre- 
sponding to the best fit model. We found the Bayes factor rises 
to 17.2. Hence, according to this corrected value for the Bayes 
factor, the ellipsoidal model is 17 times more probable than the 
spherical model meaning that the data show positive evidence 
for the deformation model. 

Furthermore, in our case we do not need to compare two 
hypothesis but we need to access the detectability of a measure- 
ment and hence we should use parameter inference instead of 
model comparison. Hence, instead of answering the question of 
whether the planet is deformed we answer the question of how 
much the planet is deformed. This latter question is best assessed 
by the analysis of the posterior probability distribution of hy 
which measures the deformation rather than by model compari- 
son. Since we found that the distribution of hy does not include 
the spherical model (hy = 0) at 30, we conclude that the de- 
formation was detected. Using the posterior distribution of hy, 
we can compute more accurate limits on the detection given that 
the distribution is not completely Gaussian. We found that hps 
is higher than 0.18 at the 99.7 confidence limit (30) and hy is 
higher than 0.03 at the 99.95 confidence limit (3.5 c). Hence, 
the detection is slightly higher than 3 c. 


5.2. Impact of limb darkening 


The model of the limb darkening can affect the measurement of 


the Love number (Akinsanmi et al.|2019| Hellard et al.|2019). 


Despite several studies on the best way to model the LD in exo- 


planet light curves (Csizmadia et a. [2013 2011), con- 


sensus still eludes us as it appears that the best model might de- 
pend on the quality of the data being analysed (e.g. 
Jordán|2015). Of the several LD parametrisations, the non-linear 
law (Claret|2000) is usually regarded as the best description of 
the stellar intensity profile (Howarth|2011); however, when fit- 
ting the parameters in the transit light curves, the correlations be- 
tween the four parameters can lead to non-physical models. Re- 
cently, the power-2 law has been shown to be 
a good balance between a small number of parameters and being 
a good approximation of the stellar intensity profiles 
[et al.|2017). Therefore, it has gained much interest helped by 
a faster algorithm (Maxted & Gill|2019). The most commonly 
used LD law is the quadratic law er 1950) due to its rela- 
tive simplicity, fast implementation, and the existence of several 
parametrisations to minimise correlations between the two pa- 
rameters (e.g. [Kipping[2013a). In addition to the choice of the 
parametrisation, it 1s also unclear if it is best to fix the LD coef- 
ficients to theoretical values based on stellar models or directly 
fit the LD in the light curves. The best approach depends not 


only on the precision of the light curves (e.g. |Espinoza & Jordan 
2015), but also on the geometry of the system (e.g. 
2011) and on whether the star is active (Csizmadia et al.|2013). 


We assessed if the LD model affected the measurement of 
the Love number by performing several tests. We compared the 
results for three different LD laws: the quadratic law which is the 
most widely used, the non-linear 4-coefficient law considered to 
be the best model, and the power-2 law which has been shown 
to give good results despite its simplicity. We fitted the LD co- 
efficients using priors derived from the stellar models. We found 
that the results depend on the priors. In particular, if the priors 


were very large, the results are independent from the stellar mod- 
els. This results in a loss of information and loss of correlations 
between the four different colours which is not ideal since they 
relate to the same star. Therefore, to try to have LD coefficients 
that are consistent for our four instrument filters, we investigated 
which were the smallest reasonable priors for the parameters for 


used the LDTK code (Parviainen & Aigrain|2015) to fit the limb 
darkening laws mentioned above for the four filters increasing 
the intrinsic uncertainty of the models, which account for the 
uncertainty in the stellar parameters, in order to obtain modelled 
law uncertainties that encompass both the PHOENIX and the 
ATLAS stellar intensity profiles. This required increasing the in- 
trinsic model uncertainties by 5 — 40x for the quadratic and the 
power-2 law. The factor is higher for the visual filters than for the 
infrared. For the non-linear law, there is no need to increase the 
model uncertainties because the four parameters give sufficient 
flexibility for the model to encompass both sets of stellar inten- 
sity profiles. The uncertainty of the modelled LD law was de- 
rived by randomly drawing LD coefficients from a Gaussian dis- 
tribution centred on the LD value and standard deviation equal to 
its uncertainty. An example, of the fit for the HST filter is shown 
in Figure[7] For CHEOPS, we could not derive the intensity pro- 
file for the ATLAS models so we used the KEPLER filter, which 
is similar to CHEOPS, as a proxy for the uncertainties needed. 
Moreover, for Spitzer, we needed to redefine the stellar radius 
for the PHOENIX models in order to match the stellar radius of 
the ATLAS models since in this case the automatic limb defini- 


tion does not give optimal results (Parviainen & Aigrain|2015; 
Espinoza & Jordán|2015). 
From Figure [7| it is clear that the PHOENIX and the AT- 


LAS models predict different stellar intensity profiles close to 
the stellar limb. It is also clear that the power-2 law matches the 
ATLAS models better and that the non-linear law matches the 
PHOENIX models better (the latter is by construction). A com- 
parison between LD derived from transit light curves and from 
theoretical models by [Espinoza & Jordán] (2015) suggested that 
the ATLAS models might be a better match to the transit fitted 
LD coefficients. However, this conclusion might depend on sev- 
eral factors that have yet to be investigated. Therefore, we expect 
that the quadratic law will be a poor description of the true stellar 
intensity profile and that the power-2 law will be a good descrip- 
tion if the ATLAS models are closer to the true stellar intensity. 
We also expect that the non-linear LD law has enough flexibility 
to match both cases. 

The uncertainties of the LD coefficients derived with the pro- 
cedure described above were used to set the priors on LD coef- 
ficients for the transit light curve fit for the spherical and the 
ellipsoidal model. The quadratic LD coefficient priors are given 
in Pod while the non-linear LD coefficient priors are given 


in Table The priors and results for the adopted model — the 
power-2 law — were already given in Section [5] We find that de- 
spite its simplicity, the power-2 law gives results in good agree- 
ment with the more complex non-linear limb darkening law. For 
the three LD laws that we tested, we obtained consistent results 
with all of the fitted parameters agreeing within 1 ø. In Table ??, 
we give the derived Love number, the significance of the detec- 
tion, and the Bayes factors. As mentioned above, for model com- 
parison, the Bayes factors should be multiplied by the prior odds 
that are very strongly in favour of the ellipsoidal model. The sig- 
nificance of the results varies slightly, but it agrees well between 
the models supporting the robustness of our results. For both the 


Article number, page 15 of 21 


A&A proofs: manuscript no. wasp103 


1.01 
0.94 
0.8 

X 0.7 

m power 2 
0.6 : quadratic 
0.51 f? . non-linear 

| — ATLAS 
0.41 —— PHOENIX 
0.0 0.2 04 06 08 1.0 
H 


Fig. 7. Stellar intensity profiles from the PHOENIX (solid purple line) 
and ATLAS (solid black line) stellar grids for the HST WFC3.IR.G141 
filter as a function of u (u = V1 -— Z2, where z is the normalised dis- 
tance from the centre of the stellar disc). We overlaid the best fit limb 
darkening models for the power-2 (dotted red), quadratic (dotted cyan), 
and non-linear (dotted green) laws. We also plotted the range of the pa- 
rameter space allowed by the limb darkening models using the derived 
parameter uncertainties after multiplying the intrinsic theoretical model 
uncertainties provided by LDTk by 40x for the quadratic model and 
10x for the power-2 model. The intrinsic uncertainties of the modelled 
grids were not changed for the derivation of the non-linear LD parame- 
ters. 


Table 9. Comparison of the derived Love number, significance of the 
detection, and Bayes factors for the three LD laws considered. 


LD law Love number Significance Bayes factor 

Power-2 law 1.597953 3c 9.1(17*) 

Quadratic 1.372959 2.30 4.6(6.6") 

Non-linear 1.6077 3.50 16(26.9*) 
Notes. 


* Corrected value of the Bayes factor as explained in Section|5.1] 


power-2 law and the non-linear law, we obtained a detection of 
the Love number of WASP-103b at more than 3 ø and consis- 
tent with each other. It is noteworthy that although the quadratic 
law provides results compatible within 1 c with the other laws, 
it yields the smallest value and the largest uncertainties for hy. 
This supports the idea that it is the worst model of the three. 
Since the three models agree well within 1 c, we conclude that 
our treatment of the limb darkening is robust and it is not biasing 
the results. 


If we increase the uncertainties of all the priors of the LD 
coefficients, for example to 0.1, we still obtain consistent values 
for hy, despite, as expected, the detection significance being re- 
duced to ~ 2c. However, we think this overestimates the true 
uncertainty of the LD, especially in the infrared where the LD 
signature is much smaller. If we use the intrinsic uncertainties 
derived from the theoretical stellar grids that are much lower than 
the ones we derived with our method (up to 40x), we find Love 
number values that are too high, indicating that the LDs were 
biasing the retrieval of the Love number. Therefore, we find that 
the best approach is to use as much prior information from the 
theoretical stellar grids as possible, while taking the differences 
associated with different models into account (in our case AT- 
LAS and PHOENIX). 
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5.8. Future prospects of measuring the tidal deformation 


Since the Love number is an important constraint for interior 
models, we tested the possibilities of constraining it better. A 
higher significance of the result would also be desirable for a 
more robust detection which requires more high signal-to-noise 
transits of WASP-103b. We simulated more seasons of CHEOPS 
observations, assuming in each season we would observe 12 
more transits. If we could obtain four more seasons of obser- 
vations (48 transits over 4 years), we would be able to mea- 
sure the Love number of WASP-103b at 4.3 o. If we could ob- 
tain six more seasons of observations (72 transits over 6 years), 
we would be able to derive hy at 5c. This would require the 
CHEOPS mission being extended. 

The extreme high precision of JWST and the fact that the 
limb darkening signature is lower in the infrared implies that the 
best chances of significantly increasing the precision of the mea- 
surement of the Love number in the near future is to combine 
our data with a transit with JWST. Since we are interested in 
maximising the cadence and the signal-to-noise of the observa- 
tions, the best would be to use the NIRSpec Prism mode which 
would enable a precision of 62 ppm/min. We simulated a transit 
of WASP-103b with NIRSpec Prism mode assuming the Love 
number derived above. We assumed a limb darkening profile 
similar to the one of Spitzer channel 1 since it has a similar wave- 
length range as the NIRSpec Prism. This simulated JWST transit 
was combined with our data and we followed the same procedure 
as above to derive the transit parameters. We obtained a 12 c de- 
tection of the Love number of WASP-103b, hy = 1.62012. This 
would be an unprecedented constraint on the Love number of an 
exoplanet which would give us strong insights into the interior of 
these planets and their similarities and differences with the Solar 
System giants. 


5.4. Measuring the Love number from planet-planet 
interaction 


HAT-P- 13 b is the only exoplanet for which the measurement of 
the Love number was confirmed. For HAT-P-13 b, the determi- 
nation of the Love number was made by an alternative method 


through the fixed point orbital eccentricity (Batygin et al.|2009; 


2012} |Buhler et al.[2016). This method, proposed 
by|Batygin et al.|(2009), is based on dynamical effects, and thus 


accesses the potential Love number, ky, instead of the radial 
Love number, Ap, as in our case. The two Love numbers are re- 
lated to each other by hy = 1+ky (e.g.|Lambeck|1980), but while 
hy solely depends on the shape of the planet, k depends on the 
knowledge of all dynamical effects in the system that can disturb 
the precession rate. It is, therefore, a much less direct method. 

Measuring the A; directly from the signature in the light 
curve only assumes that the orbit of the planet is circular and its 
rotation synchronous (Correia[2014). These two hypotheses are 
very likely for planets near the Roche limit. In contrast, many 
more assumptions are required to measure ky because the pre- 
cession rate depends on many physical effects, namely general 
relativity, rotational flattening of both the star and the planet, 
tidal deformation of both the star and planet, and secular pertur- 
bations from all the remaining planets in the system. 

applied their method to the HAT-P-13 
system, which is the only system currently known that fulfils 
the criteria of applicability. However, some assumptions were 
necessary for the derivation of the Love number. The most im- 
portant are coplanar orbits, the eccentricity is locked in a fixed 
point, and a hierarchical two-planet system. These assumptions 
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are consistent with current observations of the HAT-P-13 sys- 
tem, although they cannot be completely confirmed. If any one 
of the assumption fails, the measurement of the Love number of 
HAT-P-13 b would be biased. Nevertheless, the estimate of the 
Love number has allowed (2009), 
(2012), and (2016) to place unprecedented con- 
strains on the core mass and on the metalliticty of the planet’s 
envelop, showing the potential of the Love number to lift degen- 
erancies of the interior structure models. (2019) 
applied the same method to WASP-18Ab deriving a Love num- 
ber of kr = 0.62105, However, in this case, the cause of the 
orbital precession is not clear. 

In conclusion, the apsidal precession method allows one to 
derive precise values for the potential fluid Love number for two 
planet systems with special orbital configuration. However, the 
several assumptions of the model can have an impact on the ac- 
curacy of the measurements of the Love number. 


6. Conclusions 


We obtained 12 new high-precision transit observations of 
WASP-103b with the CHEOPS satellite to study the tidal inter- 
action with its host star. The CHEOPS data were analysed with 
a multi-dimensional GP constrained by several instrumental pa- 
rameters to correct the systematic effects due to the rotation of 
the field. We find that the roll angle, which measures the rotation 
of the field, is the instrumental parameter with higher correla- 
tion with the systematic effects. We also found that detrending 
the data with only the roll angle gives a good correction of the 
systematic noise. However, in most cases, including other in- 
strumental parameters is a better model of the systematic noise 
according to Bayesian model comparison. 

We find that a linear ephemeris is the preferred model for 
the orbital evolution of WASP-103b. However, there is a hint 
of an orbital period increase, contrary to what was expected if 
tidal decay was dominating the orbital period evolution of this 
planet. We explored scenarios that could explain a positive pe- 
riod derivative in case it is confirmed by future observations. 

One possibility is RV acceleration due to a bound compan- 
ion. If the known visual companion of WASP-103 is bound, 
it could affect the transit times of WASP-103b. To check this, 
we obtained further RV observations with CORALIE and lucky 
imaging observations with the AstroLux camera. The RV obser- 
vations are compatible with both a bound companion and a non- 
bound companion. We find an RV offset of 14 + 45 m/s between 
the previous observations and the new observations spanning 8 
years. This measured RV offset includes an unknown instrumen- 
tal offset of 14-24 m/s and the hypothetical contribution from a 
bound companion. The value of the RV offset does not exclude 
that the visual companion is bound since the RV variations over 
the 8 year timescale of the observations are expected to be less 
than 23.2 + 3.3 m/s. Although the RV required to cause the ob- 
served transit timing variations by the change in the light travel 
time is much higher (342 + 146 m/s), its high uncertainty also 
does not allow for the exclusion of this possibility. 

The new lucky imaging observations do not find the stel- 
lar companion despite the high sensitivity at the position it was 
observed before. To avoid detection, the companion star had to 
move in the direction of WASP-103 by 77 au, which is too large 
for a bound object. Therefore, the new lucky imaging observa- 
tions support the idea that the visual companion is not bound 
unless unknown systematics have affected our results. Hence, 
our data support a non-bound companion, but we recommend 


further observations to confirm these results. Long-term moni- 
toring of the RVs, as well as the new data release from Gaia, 
can help constrain the visual companion. Furthermore, high res- 
olution imaging would allow confirmation of the position of the 
visual companion and its unbound nature. 

Other possibilities to explain a positive period derivative are 
the Applegate effect and apsidial precession. However, a simpler 
explanation would be statistical artefacts. Several systematic ef- 
fects have been shown to affect the measurement from transit 
times in exoplanets (e.g. [2020). Hence, sta- 
tistical artefacts could cause the measured period to appear to be 
increasing. This is supported by the Bayesian evidence and a de- 
crease in the measured period derivative P = 3.6 + 1.6 x 107!° 
relative to previous observations ( P = 8.4 + 4.0 x 1071? - |Patra] 
et al./2020). If we assume tidal decay is dominating the period 
evolution of WASP-103b, we can place a lower limit on the tidal 
quality factor Q^ of 3.3 x 10° at 3 c, corresponding to a 99.7% 
confidence interval. This is similar to the recent limit on Q' ob- 
tained for WASP-18b 3.9x10° at a 9596 confidence interval (Ma-| 
[ciejewski et al.[2020). For these systems, longer monitoring of 
the transit times will be required in order to constrain the stellar 
tidal quality factor. 

We combined our new 12 CHEOPS light curves with pre- 
vious transit light curves obtained by HST and Spitzer to mea- 
sure the deformation of WASP-103b via its Love number. We 
re-reduced the light curves obtained with HST and Spitzer to 
correct for systematic effects since corrected light curves were 
not available in the literature. We measured the tidal deforma- 
tion of the planet directly from its distortion of the transit light 
curve. We estimated the Love number of WASP-103b to be 
hy = 1.5902, which is the first 3c detection of an exoplanet 
Love number measured directly from its deformed transit shape. 

The Love number of WASP-103b is similar to Jupiter's and 
slightly larger than Saturn's (hy = 1.39 + 0.024 — 
[2017). For a given planetary mass and radius, higher Love num- 
bers imply a metal enrichment of the envelope and hence a de- 
crease in the core mass. Our measurement of the Love number 
can be used to remove degeneracies in planetary internal models, 
allowing one to calculate the core size and the composition of 
WASP-103b (Baumeister et al./2020). Uncertainties in the limb 
darkening can influence the measurement of the Love number 
and hence we have performed a careful treatment of the limb 
darkening and several tests that indicate that our results are ro- 
bust. 

Future observations with the JWST can help to better con- 
strain the Love number of WASP-103b and gain an unprece- 
dented view of the interior of this hot Jupiter. This could help 
us to better understand these extreme systems. 
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Appendix A: Figures of the lucky imaging Table B.1. Priors for LD coefficients for the quadratic law. 
observations 


Parameter Prior 
LD1 CHEOPS  N(0.5269, 0.0218) 
LD2 CHEOPS = N(0.1279, 0.046) 
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— LD2 SP2 (0.0976, 0.0214) 
F Notes. 
= 0.00 N (a; b) is a normal distribution with mean a and standard deviation b. 
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—0.25 Table B.2. Priors for LD coefficients for the non-linear law. 
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Fig. A.2. Sensitivity curve for the AstraLux image of WASP-103. The 
1% contamination level is marked by the horizontal dotted line and the 
maximum magnitude of a blended binary to be able to mimic the transit 
of WASP-103b is marked by the horizontal dashed line. The location of 


the previously detected companion by [Ngo et al.](2016) is marked as a 


star-like symbol. 


Appendix B: Priors for the limb darkening 
coefficients 
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